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. . . We derive gravitational waveforms needed to compute the 14tli post-Newtonian (14PN) 

^SJ ' order energy flux for a particle in circular orbit around a Schwarzschild black hole, i.e. v'^^ 

beyond the leading Newtonian approximation where v is the orbital velocity of a test particle. 
, We investigate the convergence of the energy flux in the PN expansion and suggest a fitting 

■ formula which can be used to extract unknown higher order PN coefficients from accurate 

numerical data for more general orbits around a Kerr black hole. The phase difference 
between the 14PN waveforms and numerical waveforms after two years inspiral is shown to 
be about 10"^ for fi/M = 10"* and 10"^ for fi/M = 10"^ where is the mass of a compact 
object and M the mass of the central supermassive black hole. In first order black hole 
' perturbation theory, for extreme mass ratio inspirals which are one of the main targets of 

I Laser Interferometer Space Antenna, the 14PN expressions will lead to the data analysis 

accuracies comparable to the ones resulting from high precision numerical waveforms. 
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O ' Introduction. Gravitational waves emitted from a stellar mass compact object 

of mass 1 — IOOM0) orbiting a supermassive black hole of mass M(~ 10^ — 
IO^Mq) at the centers of galaxies are one of the main astrophysical sources for 
the Laser Interferometer Space Antenna (LISA). The expected number of events of 
these extreme mass ratio inspirals (EMRIs) detected by LISA during its mission 
^ I will be 10 — lOOCP (see also a recent work on the capture rate by the 2.5PN A^- 

' body simulation^). Observation of gravitational waves from EMRIs will provide 

information of the source, such as the masses, the spin of the black hole and the 
1/-^ I distribution of compact objects in the centers of galaxies. To extract such information 

^ ■ by data analysis, the accumulated difference of phase between templates and the true 

. signal during the observation should be less than one radian. The total number of 

wave cycles detected by LISA is 10^ — 10^ since LISA has maximum sensitivity around 
10~^Hz and its mission time is a few years. Thus, for LISA parameter estimation of 
EMRIs, the accuracy of theoretical waveforms should be better than 10~^. 
^ I Since the mass ratio is extreme, EMRIs can be described by black hole perturba- 

^ ■ tion theory, which uses the mass ratio /x/M as an expansion parameter (see Ref . [3]) 

■ ■ ■ for fully nonlinear numerical simulations of n/M = 1/100). At the lowest order, the 

small body traces a geodesic in the black hole geometry. Over time scales of order 
~ M'^/h (we use units G = c = 1) however, the orbit evolves adiabatically due to 
the small body's interaction with its own gravitational field, i.e. gravitational self- 
force This means that the gravitational self-force has to be taken into account 
for LISA parameter estimation of EMRIs. Second order black hole perturbation the- 
ory may also be required for the parameter estimation and is of great interest when 
the mass ratio is large. 
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Although black hole perturbation theory is powerful enough that one can com- 
pute gravitational waves accurately even in the strong field, the computational cost 
is very high to perform calculations which cover all the parameter space of EM- 
RIsP This motivates us to use analytic modeling of gravitational waves by the 
post-Newtonian (PN) theory, which uses the orbital velocity of a compact object 
V = y^M/vQ, where rg is the orbital radius, as an expansion parameter. 

In this letter, we investigate the extent of which PN order waveforms improve 
the solution for LISA parameter estimation of EMRIs by using the first order black 
hole perturbation theory. Currently, the highest available PN order for gravitational 
waveforms for EMRIs is 5.5PN {v^^ beyond Newtonian order) for a test particle in a 
circular orbit around a Schwarzschild black hole.'^' However, the convergence of the 
PN expansion becomes worse when v is larger and the phase difference between the 
PN and numerical waveforms exceeds one radian after two years inspiral.'^''^ In this 
work, we improve on the accuracy of the energy fiux to 14PN and exhibit clearly its 
closeness to a high precision numerical computation of the energy flux. We would 
like to point out that such high PN order computations have not been performed 
up till now, since the number of terms necessary to derive the PN waveforms grows 
exponentially when the PN order becomes higher. For instance, our current code 
uses 70, 3.3 x 10^, 9.0 x 10^, 1.9 x 10^, 5.6 x 10^ and 1.1 x lO^MBytes memory, 
taking seconds to half an hour, to compute multipolar waveforms for i = m = 2 
mode at 6PN, lOPN, 12PN, 14PN, 16PN and 18PN respectively. Thus, it wiU 
be difficult to obtain 19PN or higher order expressions with reasonable time on a 
personal computer by using our current code. We also suggest the fitting formula 
of the energy fiux which can be used for more general orbits around a Kerr black 
hole. One of the important consequence of this work for LISA is our demonstration 
that the phase difference after two years inspiral between our new 14PN waveforms 
and high precision numerical waveforms is negligible. This indicates that the 14PN 
expressions will lead to accuracies in LISA data analysis for EMRIs over two years 
comparable to accuracies resulting from high accuracy numerical waveforms. 

The method. The fundamental equation of black hole perturbation is the Teukol- 
sky equation, which describes the gravitational perturbation of a black hole in terms 
of the Weyl scalar '1/^ is related to the gravitational wave polarizations at 

infinity as — ihx)/2 where dot, ', denotes time derivative, d/dt. The 

Teukolsky equation can be separated if we expand in the Fourier domain using 
the —2 spin-weighted spheroidal harmonics. '3' For the case of Schwarzschild black 
hole, the radial Teukolsky equation reads 



dr \ A dr 



RimUi^) = Tem.uj{r), (1) 



where A = r{r — 2M), U{r) is the potential and T^mw is the source term depending 
on the particle's orbit. 

The radial Teukolsky equation can be solved by the Green function method. In 
this work, we use a formalism developed by Mano, Suzuki and Takasugi (MST) to 
obtain homogeneous solutions of the Teukolsky equation.® In this formalism, the 
homogeneous solutions of the Teukolsky equation are expressed in a series of hy- 
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pergeometric functions and Coulomb wave functions, which converge at the horizon 
and infinity respectively. The formalism is very powerful for the PN expansion of the 
Teukolsky equation since the series expansion is closely related to the low frequency 
expansion. The expansion coefficients of the two series aj^ are the same and satisfy 
the three-term recurrence relation 

«+l + Pn^n + 7n<-l = 0, (2) 

where u = i + 0{uj'^), < = 0{u;), = 0(1) and = 0{u)). We note that the 
parameter z/ does not exist in the original Teukolsky equation and is determined so 
that the series converges and represents a solution of the Teukolsky equation. One 
can derive the low frequency expansion of a'^ by solving Eq. ([2]) iteratively. Thus, 
if we derive aj^ up to a required order, we automatically obtain the PN expansion 
formulas up to the required order. See Ref. [T0|) for details of the formalism. 

Using the MST formalism, the energy flux absorbed into the horizon was derived 
up to relative 4PN (i.e. 6.5PN beyond Newtonian order) for a test particle in a 
circular and equatorial orbit around a Kerr black hole.'^ Gravitational wave flux to 
infinity was also computed up to 2.5PN for a test particle in slightly eccentric and 
inclined orbits around a Kerr black hole.'^''^ More recently, the MST formalism 
was applied to obtain the 5.5PN waveforms for a test particle in a circular orbit 
around a Schwarzschild black hol^ and the 4PN waveforms for a test particle in 
a circular and equatorial orbit around a Kerr black hole,'^ confirming the 5.5PN 
energy fiux in Ref. I15|) and the 4PN energy fiux in Ref. I16|) respectively. 

Results and Analysis. Once we compute the homogeneous solutions of the Teukol- 
sky equation, using the Green function method we build the solution of the Teukolsky 
equation Eq. ([T]), which is purely outgoing at infinity and ingoing at the horizon. For 
the case of a test particle in a circular orbit around a Schwarzschild black hole, the 
gravitational wave flux to inflnity is given by 

dE _ v-^ \Zemuj\'^ ,„^ 
dt ~ ^ ^ 47ra;2 ' ^"^^ 

1=2 m=-£ 

where a; = mi? is the frequency of gravitational waves, f2 = \J M/r^ is the angular 
frequency of the orbit and .^^mw is derived from the asymptotic behavior of the 
solution of the Teukolsky equation at inflnity. 

We compute gravitational waveforms which are necessary for computing the 
14PN energy flux to inflnity for a test particle in a circular orbit around a Schwarzschild 
black hole. We show only the next new 6PN terms of the energy flux, dE^^'^^/dt, 
because of space restrictions. 

dE^'^'^^ _ 2067586193789233570693 _ 246137536815857 _ 271272899815409 

dt ~ 602387400044430000 157329572400 157329572400 ^ 

246137536815857 , , , 27392 , , , 54784 , , , o 27392 2 , . ^ 

In (v) C (3) In (2) vr^ vr^ In (v) 

157329572400 ^ ^ 105 ^ ^ ^ 315 ^ ' 315 ^ ^ 

5861888 , , , , 3803225263 o 256 . 5861888 , ,,0 

+ In 2 In (v) + vr^ vr^ + (In (2 T 

11025 ^ ^ ^ ^ 10478160 45 11025 ^ ^ " 
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1465472 , ,,2 27392 o 2930944 , , , 5861888 , , , 

(In (v)) 7r^7 H In (v) 7 H In (2) 7 

11025 ^ ^ 315 ' 11025 ^ ^ ' 11025 ^ ' ' 

1465472 2 437114506833 , , , 37744140625 , , , 

7 In (3) In (5) , 

11025 ' 789268480 ^ ^ 260941824 ^ ^' 



(4) 

where C(^) is the Zeta function. We note that (ln{v))'^ term appears at 6PN. We also 
note that (ln(t)))^ and (ln(t;))^ terms appear from 9PN and 12PN respectively. One 
of the reasons is that (In(u))"' terms are produced by the FN expansion of z'^ in the 
homogeneous solution of the Teukolsky equation, where z = cor and u = £ + 0{v^). 
These (ln(f))" terms agree with Eq. (44) in Ref. I17p . which predicts (ln(?;))" terms at 
3?i-PN using the renormalization group equations. Our present work suggests that 
to obtain unknown PN terms by fits to numerical results one must include (ln(i)))'" 
terms starting from 3n-PN with the structure 

fc=0 p=0 

where [• • • ] is the floor function. Note that one can use the fitting formula Eq. ([5]) 
for more general orbits around a Kerr black hole, which are known at most up to 
4PNpJ,[i3 

We use the factorized multipolar waveforms introduced in Ref. [T8|) to compute 
the gravitational waveforms. The multipolar waveforms are decomposed into five 
factors as 

him = C''"^ Si^'^ Tim e'''- {pimY , (6) 

where /i^^''^*'^ is the Newtonian contribution to waveforms and tp denotes the parity 
of the multipolar waveforms. In the case of circular orbits, Cp = {cp = 1) when 
i + m is even (odd). The ^-th root of the amplitude pim improves the convergence 
of the PN expansion by dealing with the linear term of i in the IPN terms of the 
amplitudes. See Ref. [T5|) for the definitions of the other terms. 

The factorized multipolar waveform is the simplest and the most efficient re- 
summation technique so far. The energy flux computed from the factorized re- 
summed waveforms agrees better with the numerical results than the ones from 
Taylor-expanded waveforms and other resummed waveforms as Pade approxima- 
tion.liS 

In this letter, we show only the next new 6PN terms of /02,2) which was earlier 
computed up to 5.5PN.'I^''^ The complete expressions will be shown elsewhere.'^ 

(12) _ 313425353036319023287 6848 91592 2 

^2-2 1132319812111488000 lOS" ^ ~ 11025 ^ 

241777319107 , , , , 91592 , , , , ,,2 
- 3208936500 ^^^^^^"^ + TTo^ ^^^^^^^"^ ' 

where is 0{v"') coefficient of P2,2) eulerlog(2, v) = 7 + ln(47;) and 7 is the Euler 
constant. Again, we note that (ln(r;))2 term appears at 6PN. 
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Fig. 1. Absolute values of the difference of energy flux to infinity between numerical results and 
Taylor-expanded (left) or factorized (right) PN approximation as a function of orbital velocity. 
Note that 14PN flux converges very well even around ISCO, v = 1/VG = 0.40825. 



Fig. [T]shows the comparison of energy flux to infinity between the PN approxima- 
tions and the numerical calculation. To demonstrate the efficiency of the factorized 
waveforms, we show the results of the PN energy flux using the Taylor-expanded 
waveforms (left panel) and the factorized waveforms (right panel) in Fig. [TJ The 
numerical energy flux is obtained using the high precision code in Ref. flUti . The 
accuracy of the numerical calculation is better than 10"^*^ if we set the maximum 
value of i to 20 in Eq. (jS]). The n-PN flux needs £ up to n + 2. The agreement of the 
total energy flux between the PN and the numerical results becomes better when 
the PN order is higher even around the innermost stable circular orbit (ISCO) (see 
Ref.EU for the calculation beyond ISCO). 

To investigate quantitatively the PN order needed for LISA parameter estima- 
tion of EMRI, we compare the phase difference during two years quasi-circular inspi- 
ral between the factorized PN waveforms and the numerical waveforms. Following 
Ref. [7]), we examine two systems, which are named System-I and System-II. System-I 
has masses (M,/i) = (1O^,1O)M0 and insphals from vq ~ 29.34M to tq ~ 16. IM, 
whose frequency sweeps from /qw — 4 x 10~^IIz to /gw — 10~^Hz. System-II has 
masses (M, /i) = (1O^,1O)M0 and starts inspiral from tq ~ 10. 6M to tq ~ 6.0M, 
whose frequency sweeps from /qw — 1-8 x lO^'^Hz to /gw — 4.4 x 10~'^Hz. System-I 
(II) has ~ 1 X 10^ (~ 5 X 10^) rads of orbital phase during its inspiral and corresponds 
to the early (late) inspiral phase of an EMRI in the frequency band of LISA. 

Fig. [2] shows the comparison of the phase between the dominant mode /i2,2 of the 
factorized PN and the numerical calculation (see Ref. [22]) for the calculation of the 
phase). For System-I (II), the absolute values of the phase difference between the 
factorized PN waveforms and the numerical waveforms after two years inspiral are 
about 3.6 (4.0 x 10^), 0.38 (1.2 x 10^), 5.4 x 10"^ (10), 3.4 x 10"^ (0.57), 1.8 x lO"*' 
(0.11) and 2.5 x 10"^ (8.9 x 10"^) rads for 5.5PN, 6PN, 8PN, lOPN, 12PN and 14PN 
respectively (The relative error of the amplitude between the factorized 14PN and 
the numerical waveforms is 6.6 x 10~^^ (7.1 x 10"'') for System-I (II)). 

The strongest EMRI events detected by LISA may have the signal to noise ratio 
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Fig. 2. Absolute values of the phase difference due to two years inspiral between the factorized 
PN and the numerical waveforms for ^2,2 as a function of time in month. The left panel 
shows the dephase for {M,fi) = (10'^,10)Mq, which starts from ro ~ 29.34M to ro ~ 16. IM 
and sweeps frequencies /gw € [4 x 10"'', 10~^]Hz. The right panel shows the dephase for 
{M,n) = (1O^,1O)M0 from ro — 10. 6A/ to ro — 6.0AI with associated frequencies /gw G 
[1.8 X 10~^,4.4 X 10~^]Hz. The left (right) panel represents the early (late) inspiral phase in 
the LISA band. If the dephase is less than 10 milliradians, the PN waveforms will provide the 
parameter estimation of EMRI comparable to the one resulting from numerical waveforms. 



up to p 100.'^ Thus, LISA can measure phase difference of the order of 1/p ~ 10 
mihiradians by matched filtering.'S' This suggests that the 8PN (14PN) waveforms 
win lead to LISA parameter estimation of System-I (II) comparable to the one using 
numerical waveforms. Since System-II represents the inspiral in the most strong-field 
of a Schwarzschild black hole, the 14PN waveforms will provide LISA parameter 
estimation of EMRI comparable to the one resulting from numerical waveforms. 

Conclusion and Discussion. Using the first order black hole perturbation theory, 
we have computed gravitational waveforms consistent with the 14PN energy flux for 
a test particle in a circular orbit around a Schwarzschild black hole. The high PN or- 
der computation has been performed systematically using the MST formalism.lS We 
provide a fitting formula of the energy flux, Eq. ([5]), for more general orbits around 
a Kerr black hole. Comparing the energy flux with a high precision numerical com- 
putation, we investigated the approach of the PN expansion towards the numerical 
results and found that the PN expansion converges well even at ISCO. The phase 
difference between our new analytic waveforms and the numerical waveforms after 
two years of inspiral is so negligibly small that it does not lead to any discrepancy 
in LISA data analysis. Thus, using first order black hole perturbation theory, one 
can derive sufficiently high PN order waveforms analytically and build the template 
banks to span the parameter space efficiently without recourse to a more expensive 
numerical computation. 

The EMRIs contain more generic systems such as eccentric and inclined orbits of 
a test particle around a Kerr black hole (see Ref. for an application to a spinning 
particle). The orbital velocity at ISCO can be larger when the black hole is rotating. 
If the orbits become more generic, one needs a larger number of wave modes, which 
will contain higher frequency contributions than the case in this work. Thus, for 
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more generic orbits around a Kerr black hole we have to compute higher PN orders 
than 14PN in order to achieve good agreement with the numerical waveforms. Since 
the formalism adopted in this work is systematic and does not have theoretical issues 
to compute higher PN order, in principle one can compute sufficiently high PN order 
waveforms even for generic orbits. It will be possible to perform sufficiently high PN 
order computation for the case of a test particle in slightly eccentric and inclined 
orbits around a Kerr black hole. However, it may be difficult to compute them for 
the case of a test particle in large eccentric and inclined orbits around a Kerr black 
hole in reasonable time. Moreover, conservative effects of self-force, which are not 
included in this work, may contribute ~ 20 radians to the phase of waveforms during 
inspiraPS' (however see Table II in Ref . I26p , which suggests that conservative effects 
up to 9PN is sufficient to reduce the phase error to 10 milliradians by investigat- 
ing 18PN scalar self- force). Another approach to discuss generic orbits may be the 
effective-one-body formalism which can determine unknown terms in PN approx- 
imation and self-force effects by calibrating them with numerical calculation.'^''^ 
The present work has implications for this approach too. For the calibration, we 
recommend the use of our proposed fitting formula Eq. dSJ dealing with (In(u))"' 

terms appearing from 3n-PN. 
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